Potential Vorticity Formulation of Compressible Magnetohydrodynamics 
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Compressible ideal magnetohydrodynamics (MHD) is formulated in terms of the time evolution 
of potential vorticity and magnetic flux per unit mass using a compact Lie bracket notation. It is 
demonstrated that this simplifies analytic solution in at least one very important situation relevant 
to magnetic fusion experiments. Potentially important implications for analytic and numerical 
modelling of both laboratory and astrophysical plasmas are also discussed. 

PACS numbers: 52.30.Cv, 52.55.Fa, 96.60.Q- 



INTRODUCTION 

Ideal MHD is a model for magnetised plasma where 
the coUisionality is low, so that dissipative effects can 
be neglected, yet where the charged particles still inter- 
act sufficiently strongly via the electromagnetic field they 
can be treated as a single fluid. The ideal MHD model is 
applied to a wide range of laboratory and astrophysical 
situations, where there are long periods of relative qui- 
escence in which Maxwellian particle distributions can 
be approached, interrupted by often violent transients. 
Ideal MHD instabilities are thought to be implicated in 
the triggering of the sawtooth crash phenomenon in toka- 
mak magnetic fusion experiments and flaringin the solar 
and stellar context, see textbooks such as The for- 
mer is important as it limits the performance of devices 
ultimately intended to generate nuclear power, and the 
latter is implicated in the generation of solar magnetic 
storms which can disrupt terrestrial power grids, naviga- 
tion and communication systems. Both these topics are 
presently the subject of intensive investigation, magnetic 
fusion as the multi-billion dollar ITER tokamak enters 
the construction phase, whereas multiple satellite mis- 
sions are collecting data on solar and stellar magnetic 
fields. 

It is often mathematically convenient when employing 
ideal MHD, to assume that the plasma fluid is incom- 
pressible, but the reality in the above-mentioned situa- 
tions is that the plasma density varies by one or more or- 
ders of magnitude over the region of interest. This work 
presents what is believed to be a novel, mathematically 
convenient formulation of compressible MHD. 

The equations of ideal MHD as usually formulated are 
well-known and are to be found in many textbooks, see 
eg- [H, §4.3]. As explained there, the problem admits a 
variational formulation which is of great utility for prac- 
tical stability analysis, and a functional Hamiltonian for- 
mulation in terms of Lie derivatives [2] , of great theoreti- 
cal importance for understanding stability and evolution. 
More direct approaches to ideal MHD stability are also 
now used [E § 6] , and the results presently to be described 
are more relevant to the latter school. 

The potential vorticity is the ordinary vorticity lj of 



the plasma (the curl of the mean flow U of ions and elec- 
trons), divided by the mass density p, ie. lj — u;/p. The 
possibility of combining the equation for the time evolu- 
tion of vorticity with that for density evolution to give a 
simple equation for the rate of change of potential vor- 
ticity, was first realised for a classical fluid by Hclmholtz 
as described by [3|, § 146] in the mid-19th Century. In 
the mid-20th Century, Walcn, according to [1, § 4-2] was 
the first to realise that a mathematically identical rela- 
tion governed the evolution of the magnetic flux per unit 
mass B = B//9 where B is the magnetic fleld. For incom- 
pressible plasma, Arnold & Khesin §I.10.C] combined 
these results in late-20th Century to give an elegant for- 
mulation of ideal MHD in terms of Lie brackets of vector 
flelds. The Lie bracket is here the generalisation to arbi- 
trary vector flelds of the 'flux-freezing' operator, ie. the 
operator which determines the advection of divergence- 
free (solenoidal) flelds B and a; 0, §3.8]. The novelty of 
the present work is to extend this formalism to compress- 
ible MHD and explore the implications. In particular, the 
peculiar, coordinate invariant nature of the Lie bracket 
makes it easy to generalise solutions to arbitrary geome- 
try in some cases, both analytically and numerically. 

The next section contains a detailed mathematical 
derivation of the key formula. A discussion of the im- 
plications for analytic and numerical solution follows, 
and flnally some important possible applications are sum- 
marised. 



MATHEMATICS 

In terms of the operators of Classical Vector Mechan- 
ics, the Lie derivative of a vector can be deflned as: 



£„(v) = Vx(uxv)-uV-v + vV 



(1) 



which will help explain the equivalence with the vector 
advection operator, the flrst term on the right. Indeed, 
Walen's result for magnetic induction in a perfectly con- 
ducting medium is 



dB 

'dt 



/:u(B) 



(2) 
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Introducing component notation for vectors in gen- 
eral non-orthogonal coordinate systems, as described in 
many textbooks e.g. 0, it turns out that the Jacobians 
thereby introduced (of the co-ordinate transformation 
from Cartesians), cancel among the terms in Eq. ([T|), so 
that 



(3) 



where , are the contravariant components of the 3- 
vectors u, v respectively, and the summation convention 
is implied. It follows that 



'Cu(v) = -£v(u) ^ -[u,v] 



(4) 



where [•, •] denotes the Lie bracket of Schutz [81]. 

It will be now be proved that the equation for the evo- 
lution of potential vorticity in compressible ideal MHD 
may be written 



du3 

'dt 



(5) 



where the potential current J = V x B/p. The customary 
vorticity equation in ideal MHD is 

9w „ , VpxVp „ /JxB\ 
_ = V X (U X + ^ . + V X (6) 



dt 



where vorticity a; = pu) = V x U, and current 3 = p3 = 
V X B = V X (pB). When proceeding further, it is con- 
venient and often physically justifiable, by a barotropic 
or isentropic assumption, to neglect the term in the pres- 
sure p, and if not, the resulting additional term is easily 
representable in general geometry. 

It follows that to establish the equivalence of 
Eqs dU and (|6]), it is necessary to show that A = 0, 
where 



, 1^ /BxJ 

A = -V X 

P VP 



(7) 



Now, Eq. ([7]) is a vector equation, so validity in any co- 
ordinate frame implies validity in all, hence it is sufficient 
to establish the result in Cartesian coordinates, where 

1, 



A = -V X (pB X J) +B • VJ - J • VB (8) 



The curl term may be expanded using the identity 

-V X (pv) = Rxv-|-Vxv 
P 



(9) 



where R = ^p/p- Setting v = B x J, and expanding the 
resulting curl-cross operation, there is cancellation of the 
two terms from the Lie derivative, leaving 



A = BV • J - JV • B + R X (B X J) 



(10) 



Since V • J = 0, it follows that 



V J = R J 



and likewise since V • B = 0, 



V B = -R B 



(11) 



(12) 



Substituting Eq. (HII) and Eq. (IT2|) in Eq. (HO]), and ex- 
panding the last term as dot products, shows that, as 
required A = 0. 

The set of evolution equations is completed by mass 
conservation 



dp 

dt 



-V • (pU) 



(13) 



This does not involve a vector Lie derivative, but, using 
the standard expression for the divergence operator in 
general curvilinear coordinates, it may be written 



dp 

dt 



1 djpy^W^) 

y/g dx^ 



(14) 



where ^ is the Jacobian and the gik is the metric tensor, 
which upon introducing p — p^ may be written 



dp 

dt 



d(pU^ 
dx^ 



(15) 



provided that ^ docs not change with time. Like the 
neglect of the pressure term above, this latter inessential 
assumption is often physically reasonable. 

Unfortunately, the ideal MHD equations are here com- 
pleted by the two definitions of potential vorticity and 
potential current, which do explicitly contain metric in- 
formation, viz. 



pw* 



dx'^ 



and 



pr^e^'^^i^pB 



dx^" \^ 



gin 



(16) 



(17) 



In the above, e = eiu is the alternating symbol, taking 
values 1, —1 or 0, depending whether {ikl) is an even, 
odd or non-permutation of (123). Finally, note that Eq. 
^ and Eq. (fT5|) together ensure that V • B = 0, only if 
initially 



a(pB^^ 
dx^ 



= 



(18) 



SOLVING THE NEW SYSTEM 



The new model system for ideal barotropic compress- 
ible MHD evolution consists of Eq. (O, Eq. (g)), Eq. p^ . 
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Eq. ([T6)) and Eq. (IT7)) . The simplification of tlie first three 
has been gained at the expense of complicating the last 
two 'static' relations. Nonetheless, evolution equations 
are harder to treat numerically, because any errors in the 
discretisation tend to combine over time. Moreover, it 
will be evident that problems solved in Cartesian geom- 
etry will test all aspects of the coding of the evolution- 
ary equations. Thus, there is considerable computational 
advantage to be gained. There is obviously the concern 
that the magnetic field computed may not be accurately 
solenoidal, but this is an issue for many other discretisa- 
tions also. The main difficulty is in the inversion of Eq. 
(fT6|) to give the velocity field U corresponding to a freshly 
evolved potential vorticity (since B itself is evolved, Eq. 
([T7| does not need to be inverted). However, this inver- 
sion, together with the computation of the irrotational 
part of U, is a classical hydro dynamical problem, and a 
variety of strategies may be found in the literature. On 
present machine architectures, introducing the vector po- 
tential for velocity then solving the coupled system Eq. 

and Eq. by a pseudo-timestepping algorithm 
is probably to be preferred. Similar numerical solution 
strategies were successfully employed in electromagnetics 
by the current author and collaborators ■ Vorticity 

formulations are common in plasma modelling as they 
are helpful in several physically relevant limits, and in 
particular, a vorticity formulation has been used success- 
fully in nonlinear, compressible MHD [V\\. 

Turning to analytic results, first consider MHD equi- 
librium solutions with no time dependence and U = 0, 
implying £3 (J) = 0. In the case of force-free fields, 
meaning J oc B, substituting J = AB in the Lie deriva- 
tive in component form show this is a solution provided 
B • VA = 0, i.e. exactly the same constraint on A that 
follows from the solenoidal constraint on B and J when 
seeking the solution J = AB. Hydromagnetic force- free 
solutions, with the additional constraint that U = A2B, 
now cease to exist however, because U is not solenoidal 
unless the flow is incompressible. 

Moving now to time dependent solutions, interest at- 
taches to the 'flux compression' solution § 4.6], which 
is postulated on purely kinematic grounds (i.e. from Eq. 
and which may be written 

B = c{0,0,p{x,y,t)) (19) 

for a compressible flow U with density p provided that 
U = {Ux{x,y,t),Uy{x,y,t),Q). Here, c is an arbitrary 
constant and (x, y, t) are the usual Cartesian coordinates. 
This solution is of practical importance for fusion exper- 
iments, where external magnets are used to generate a 
time dependent flux designed so as to compress plasma 
'frozen' to it. It is easy to establish that if B = cpz then 
J X B/c^ ~ (Vp X z) X pz = -V(ip^), and so there 
are compressible MHD solutions of the form Eq. (IT5|) . for 
2-D solutions of compressible hydrodynamics compatible 




FIG. 1: An, in effect unmagnetised, compressible flow is 
shown. The motion consists of rolls swirling about a B-field 
aligned with the z-axis, with arrow-heads indicating the sense 
of motion of each eddy. 

with an additional pressure gradient of this form. One 
possibility is illustrated in simple geometry in Figure [T] 

The simple form of the new evolution equations enables 
a generalisation of the flux-compression solution to gen- 
eral curvilinear coordinates. It is important to emphasise 
that the following is not simply re-expressing B = pz in 
different coordinate systems, nor is there a loss of gener- 
ality in choosing units for density such that c = 1. The 
obvious generalisation is to take = 1 = = 0); 
implying a 2-D density to ensure a solenoidal B, since Eq. 
([T8l) requires dp/dx^ = 0. The next step is to ensure that 
£3 (J) = 0, which as may be seen using the coordinate 
form Eq. ([3|), simply requires dJ^ /dx^ = 0. Similarly 
£u(B) = may be satisfied by a flow with dU^ /dx^ = 
(note that C/'^ ^ is therefore allowed). From the 'static' 
relations, it will be seen that a solution with indepen- 
dent of x^ is possible provided dgik/dx^ — 0. Put in 
the language of differential geometry [sj, §3.11], if B is a 
Killing vector, there is a flrrx-compression solution. 

Further to explore the implications of this, introduce 
generalised toroidal coordinates {g^s^w] (cf. {r,9,4>) as 
commonly employed in plasma physics [7]) so that 

X = {x,y,z) = {Rc cos w, Rc sin w, ipc sin s) , (20) 

where 

Rc = Ro + iJc{Q,s,w) cos s (21) 

It will be seen that ipcio^ s^w) = const, as g varies form 
a set of nested toroidal surfaces with major axis i?o- 
Introduce helical coordinates (u, v) on each surface, so 
that s = u — v/q{ip), w — V + u/q{ip), and write 
i){Q,u,v) = iI}c{q, s,w). Suppose that -0 is rotation- 
ally symmetric about the z-axis and satisfies the Grad- 
Shafranov equation, ie. if) is a, flux function for an equi- 
librium magnetic field, then the curves of Eq. (pn|) as 
V varies at constant u and tp are equivalent to lines of 



4 



the equilibrium field with helical twist q{ip). (Note that 
u and V need only be suitably periodic functions of the 
regular toroidal angles 9 and (j)- To define an equilibrium 
fully requires defining these functions, but this is inessen- 
tial for what follows.) The metric tensor in a coordinate 
system {x^,x'^,x^) is given by 



Taking = {g,u,v) and using suffix ,i to de- 

note differentiation with respect to x*, the components 
of gik are straightforwardly calculated as 

9ik = ip,i'>P,k + i^'^s^iS^k + (Ra + ip cos s)'^w^iW^k (23) 

and when q is constant, s^i — (0,1,-1/(7), w,i = 
(0,l/g,l). 

This x'^ coordinate system has been chosen so that the 
equilibrium field expected in the tokamak confinement 
device may be expressed as -B'^ = 1 {B^ — = 0), but 
it will be seen that in general, the metric tensor does 
depend on x^ — v through s = u — v/q. By inspection, 
however, in the limit when q is large, gik depends only on 
x^ and x^. Hence a purely toroidal field, ie. one tangent 
to circles about the major axis x = ?/ = of the torus, 
allows for fiux freezing solutions. Further, when ip/ Rq is 
small, the s-dependence of gik is weak, so a helical field 
in a torus with relatively large major radius is also in 
this category. The preceding limits illustrate two of the 
Killing vector solution symmetries [1, § 5.2.4] (the third 
is simply invariance in a Cartesian coordinate). 

Other possibilities for new analytic solutions outside 
of i?'^ = 1 are opened up when it is realised that the 
helical field considered above is just one example of the 
use of Clebsch variables 0, § 5] to represent a solenoidal 
vector field as a single contravariant component. Alter- 
natively, the vector potential may be introduced, leading 
to an interesting calculus involving R, consistent with the 
fact that exponentially varying density profiles (implying 
constant R) are often studied analytically. 

APPLICATIONS 

For magnetic fusion physics, the above, new analytic 
flux-compression solutions represent a possible nonlinear 
development of interchange modes ^3|, §12.1.2]. They 
would seem to represent an efficient and rapid means 
whereby mass (and hence heat) might escape from a dis- 
charge, hence might be implicated in situations where 
there is rapid transient cooling, such as the sawtooth 
crash in the centre of the tokamak discharge (-0 small), 
and ELMs (Edge Localised Modes) in divertor discharges 
(large q limit). The preceding section has also speculated 
that the new formalism could be used efficiently to sim- 
ulate ideal MHD evolution of discharges in generalised 



coordinates, say defined by an arbitrary MHD equilib- 
rium. 

In astrophysics, observed magnetic fields usually ex- 
hibit a significant degree of disorder, so it is unclear how 
important the new flux-compression solutions might be, 
as they rely on at least a degree of coordinate invari- 
ance. It is speculated that, in stars with a strong inter- 
nal toroidal field (such as the Sun is believed to possess), 
the rotationally symmetric solution might help model the 
convection pattern, accounting for the largely latitudinal 
variation of the solar differential rotation. Regardless, it 
should be helpful that, in the new equations, the field 
geometry appears only in the state equations. It will for 
example, be simpler to generate more realistic solutions 
from symmetric ones by varying gn^ starting with the unit 
tensor. This could be useful, say, for modelling sunspot 
penumbrae both analytically and computationally, since 
there the magnetic field is predominantly directed radi- 
ally outwards in the horizontal direction. 
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